#程序文件 - mme2024-example-7-11.py

#首先载入数值计算和函数画图的numpy 和pylab. 
import numpy as np
import pylab as plt

#然后载入数据，并将五个观测点的坐标计算成线性方程组的系数矩阵。
a=np.loadtxt('data7-11.txt')
x=a[0]
y=a[1]
A=np.vstack([x**2,x*y,y**2,x,y]).T  #线性方程组系数矩阵

#然后写出常数项的列向量，使用numpy的矩阵求逆函数 inv 与矩阵乘法函数 @, 求解线性方程组。
b=-np.ones(5)   #线性方程组的常数项
c=np.linalg.inv(A)@b  #解线性方程组拟合参数
print('coefficients are: \n',np.round(c,4))

#接下来开始画图，先将椭圆方程的左边定义成一个二元函数
#f(x,y) = c_0x^2+c_1xy+c_2y^2+c_3x+c_4y
f=lambda x,y: c[0]*x**2+c[1]*x*y+c[2]*y**2+c[3]*x+c[4]*y

#然后定义网格点，计算二元函数在每个网格点的函数值。
x=np.linspace(3,8,100)
y=np.linspace(-1,5,100)
x,y=np.meshgrid(x,y)
z=f(x,y)

#最后画出二元函数 $f(x,y)$ 取值为 $-1$ 的等值线。
#避免使用 LaTeX 渲染文本
import matplotlib
matplotlib.rcParams['text.usetex'] = False
plt.contour(x,y,z,[-1])   #画高度为-1的等高线
plt.xlabel('x')
plt.ylabel('y',rotation=1)
plt.savefig('mme2024-example-7-11.png')  #保存图像


